// do /Users/shoude/Dropbox/eegap/EEgap_data_code_heter_SM/Code/plotting/plot_Fox_contour_3dim_bs_robust.do
//Options: discrete 44000 16 1


pause on

global  pathresults = "/Users/shoude/Dropbox/eegap/EEgap_data_code_heter_SM/Results_estimation"
global  Est_pathresults = "/Users/shoude/Dropbox/eegap/EEgap_data_code_heter_SM/Results_estimation/demand_est_grID" 




local FOX_model = `"`1'"'   // discrete 
//discretev5, discretev6, discretev7, discretev8, discrete_3tau 
local subsample = `"`2'"'   // Use 44000 
local inc       = `"`3'"'   // 16
local set_seed  = `"`4'"'   // 1

insheet using $Est_pathresults/Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, clear


 forvalues i=1(1)16 { 
	capture destring weight`i', replace force 
 }
 
 collapse(sum) weight*, by(theta eta)
 
//forvalues i=1(1)10 {
//forvalues i=1(1)7 { 
forvalues i=1(1)16 { 
 capture gen eta_avg_tmp`i' = eta*weight`i' 
 capture gen theta_avg_tmp`i' = theta*weight`i' 
 capture gen weight_100`i' = round(100*weight`i')
 capture egen eta_avg`i' = sum(eta_avg_tmp`i') 
 capture egen theta_avg`i' = sum(theta_avg_tmp`i')
}  

drop eta_avg_tmp*
drop theta_avg_tmp*

reshape long weight weight_100 eta_avg theta_avg, i(theta eta) j(bs)

keep if weight_100>0
drop if weight_100==. 
sort eta theta
by eta theta: egen mean_weight=mean(weight_100)
by eta theta: egen sd_weight=sd(weight_100)
by eta theta: egen count_weight=count(weight_100)
gen  se_weight=  sd_weight/(count_weight^0.5)
gen ICl_weight=  mean_weight - 1.96* se_weight
gen ICh_weight=  mean_weight + 1.96* se_weight


preserve 
collapse(mean) mean_weight se_weight ICl_weight ICh_weight,by(eta theta)

	mvencode se_weight,mv(0) over
	replace  ICl_weight = mean_weight if ICl_weight == .
	replace  ICh_weight = mean_weight if ICh_weight == .
	gen tstats = mean_weight/se_weight
	replace tstats = 0 if tstats == .
	order eta theta mean_weight se_weight ICl_weight ICh_weight tstats 
	sort eta theta
outsheet using $pathresults/plot_jointpdf_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, comma noquote noname replace 

restore

collapse(firstnm) eta_avg theta_avg,by(bs)
	gen eta_avg_sd = eta_avg
	gen theta_avg_sd = theta_avg
	gen eta_avg_n = eta_avg
	gen theta_avg_n = theta_avg
collapse(mean) eta_avg theta_avg (sd) eta_avg_sd theta_avg_sd  (count)	eta_avg_n  theta_avg_n
	gen eta_avg_se = eta_avg_sd / (eta_avg_n^0.5)
	gen theta_avg_se = theta_avg_sd / (theta_avg_n^0.5)
	keep eta_avg theta_avg eta_avg_se theta_avg_se
outsheet using $pathresults/means_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, comma noquote noname replace 


//theta
insheet using $Est_pathresults/Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, clear

 forvalues i=1(1)16 { 
  capture destring weight`i', replace force 
 }
 
 collapse(sum) weight*, by(theta)

//forvalues i=1(1)10 {
//forvalues i=1(1)7 { 
forvalues i=1(1)16 { 
 capture gen eta_avg_tmp`i' = eta*weight`i' 
 capture gen theta_avg_tmp`i' = theta*weight`i' 
 capture gen weight_100`i' = round(100*weight`i')
 capture egen eta_avg`i' = sum(eta_avg_tmp`i') 
 capture egen theta_avg`i' = sum(theta_avg_tmp`i')
}  

//drop eta_avg_tmp*
drop theta_avg_tmp*

reshape long weight weight_100 theta_avg, i(theta) j(bs)

keep if weight_100!=0
sort theta
by theta: egen mean_weight_theta=mean(weight_100)
by theta: egen sd_weight_theta=sd(weight_100)
by theta: egen count_weight_theta=count(weight_100)
gen  se_weight_theta=sd_weight_theta/(count_weight_theta^0.5)
gen ICl_weight_theta=  mean_weight_theta-1.96* se_weight_theta
gen ICh_weight_theta=  mean_weight_theta+1.96* se_weight_theta

collapse(mean) mean_weight se_weight ICl_weight ICh_weight,by(theta)
	mvencode se_weight,mv(0) over
	replace  ICl_weight = mean_weight if ICl_weight == .
	replace  ICh_weight = mean_weight if ICh_weight == .
	gen tstats = mean_weight/se_weight
	replace tstats = 0 if tstats == .
	order theta mean_weight se_weight ICl_weight ICh_weight tstats 
	sort theta
outsheet using $pathresults/plot_marginaltheta_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, comma noquote noname replace 


//eta
insheet using $Est_pathresults/Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, clear

 forvalues i=1(1)16 { 
  capture destring weight`i', replace force 
 }
 
 collapse(sum) weight*, by(eta)

//forvalues i=1(1)10 {
//forvalues i=1(1)7 { 
forvalues i=1(1)16 { 
 capture gen eta_avg_tmp`i' = eta*weight`i' 
 capture gen theta_avg_tmp`i' = theta*weight`i' 
 capture gen weight_100`i' = round(100*weight`i')
 capture egen eta_avg`i' = sum(eta_avg_tmp`i') 
 capture egen theta_avg`i' = sum(theta_avg_tmp`i')
}  

drop eta_avg_tmp*
//drop theta_avg_tmp*

reshape long weight weight_100 eta_avg, i(eta) j(bs)

keep if weight_100!=0
sort eta
by eta: egen mean_weight_eta=mean(weight_100)
by eta: egen sd_weight_eta=sd(weight_100)
by eta: egen count_weight_eta=count(weight_100)
gen  se_weight_eta=sd_weight_eta/(count_weight_eta^0.5)
gen ICl_weight_eta=  mean_weight_eta-1.96* se_weight_eta
gen ICh_weight_eta=  mean_weight_eta+1.96* se_weight_eta

collapse(mean) mean_weight se_weight ICl_weight ICh_weight,by(eta)
	mvencode se_weight,mv(0) over
	replace  ICl_weight = mean_weight if ICl_weight == .
	replace  ICh_weight = mean_weight if ICh_weight == .
	gen tstats = mean_weight/se_weight
	replace tstats = 0 if tstats == .
	order eta mean_weight se_weight ICl_weight ICh_weight tstats 
	sort eta
outsheet using $pathresults/plot_marginaleta_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, comma noquote noname replace 




//theta: [<0, 0-0.75, 0.75-1.25, >1.25]
insheet using $Est_pathresults/Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, clear

 forvalues i=1(1)16 { 
  capture destring weight`i', replace force 
 }

 gen theta_4q =. 
 replace theta_4q = 1 if theta < 0
 replace theta_4q = 2 if theta >= 0 & theta<0.75
 replace theta_4q = 3 if theta >= 0.75 & theta < 1.25
 replace theta_4q = 4 if theta >= 1.25

 collapse(sum) weight*, by(theta_4q)

//forvalues i=1(1)10 {
//forvalues i=1(1)7 { 
forvalues i=1(1)16 { 
 capture gen eta_avg_tmp`i' = eta*weight`i' 
 capture gen theta_avg_tmp`i' = theta*weight`i' 
 capture gen weight_100`i' = round(100*weight`i')
 capture egen eta_avg`i' = sum(eta_avg_tmp`i') 
 capture egen theta_avg`i' = sum(theta_avg_tmp`i')
}  

//drop eta_avg_tmp*
drop theta_avg_tmp*

reshape long weight weight_100 theta_avg, i(theta_4q) j(bs)

keep if weight_100!=0
sort theta_4q
by theta_4q: egen mean_weight_theta=mean(weight_100)
by theta_4q: egen sd_weight_theta=sd(weight_100)
by theta_4q: egen count_weight_theta=count(weight_100)
gen se_weight_theta=sd_weight_theta/(count_weight_theta^0.5)
gen ICl_weight_theta=  mean_weight_theta-1.96* se_weight_theta
gen ICh_weight_theta=  mean_weight_theta+1.96* se_weight_theta

collapse(mean) mean_weight se_weight ICl_weight ICh_weight,by(theta_4q)
	mvencode se_weight,mv(0) over
	replace  ICl_weight = mean_weight if ICl_weight == .
	replace  ICh_weight = mean_weight if ICh_weight == .
	gen tstats = mean_weight/se_weight
	replace tstats = 0 if tstats == .
	order theta mean_weight se_weight ICl_weight ICh_weight tstats 
	sort theta
outsheet using $pathresults/plot_marginaltheta_4q_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, comma noquote noname replace 





//Alternative cut-offs for theta: based on the following discount rate: below 2% 2-12% 12%-very large, wrong sign
//theta: [<0, 0-0.75, 0.75-1.25, >1.25]
insheet using $Est_pathresults/Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, clear

 forvalues i=1(1)16 { 
  capture destring weight`i', replace force 
 }

 gen theta_4q =. 
 replace theta_4q = 1 if theta < 0
 replace theta_4q = 2 if theta >= 0 & theta < 0.6202
 replace theta_4q = 3 if theta >= 0.6202 & theta < 1.2825
 replace theta_4q = 4 if theta >= 1.2825

 collapse(sum) weight*, by(theta_4q)

//forvalues i=1(1)10 {
//forvalues i=1(1)7 { 
forvalues i=1(1)16 { 
 capture gen eta_avg_tmp`i' = eta*weight`i' 
 capture gen theta_avg_tmp`i' = theta*weight`i' 
 capture gen weight_100`i' = round(100*weight`i')
 capture egen eta_avg`i' = sum(eta_avg_tmp`i') 
 capture egen theta_avg`i' = sum(theta_avg_tmp`i')
}  

//drop eta_avg_tmp*
drop theta_avg_tmp*

reshape long weight weight_100 theta_avg, i(theta_4q) j(bs)

keep if weight_100!=0
sort theta_4q
by theta_4q: egen mean_weight_theta=mean(weight_100)
by theta_4q: egen sd_weight_theta=sd(weight_100)
by theta_4q: egen count_weight_theta=count(weight_100)
gen se_weight_theta=sd_weight_theta/(count_weight_theta^0.5)
gen ICl_weight_theta=  mean_weight_theta-1.96* se_weight_theta
gen ICh_weight_theta=  mean_weight_theta+1.96* se_weight_theta

collapse(mean) mean_weight se_weight ICl_weight ICh_weight,by(theta_4q)
	mvencode se_weight,mv(0) over
	replace  ICl_weight = mean_weight if ICl_weight == .
	replace  ICh_weight = mean_weight if ICh_weight == .
	gen tstats = mean_weight/se_weight
	replace tstats = 0 if tstats == .
	order theta mean_weight se_weight ICl_weight ICh_weight tstats 
	sort theta
outsheet using $pathresults/plot_marginaltheta_4q_2_12_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, comma noquote noname replace 



//Compute distribution for different lifetimes: 8, 10, and 12..

insheet using $Est_pathresults/Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, clear

 forvalues i=1(1)16 { 
  capture destring weight`i', replace force 
 }
 
 gen rho_5=1/(1+0.05)
 gen LLC_5_18=rho_5*(1-rho_5^18)/(1-rho_5)

 gen LLC_5_4  = rho_5*(1-rho_5^4)/(1-rho_5)
 gen LLC_5_8  = rho_5*(1-rho_5^8)/(1-rho_5)
 gen LLC_5_10 = rho_5*(1-rho_5^10)/(1-rho_5)
 gen LLC_5_12 = rho_5*(1-rho_5^12)/(1-rho_5)
 gen LLC_5_24 = rho_5*(1-rho_5^24)/(1-rho_5)

 gen theta_m_4_L    = theta*LLC_5_18/LLC_5_4
 gen theta_m_8_L    = theta*LLC_5_18/LLC_5_8
 gen theta_m_10_L   = theta*LLC_5_18/LLC_5_10
 gen theta_m_12_L   = theta*LLC_5_18/LLC_5_12
 gen theta_m_24_L   = theta*LLC_5_18/LLC_5_24
 
 gen rho_2=1/(1+0.02)
 gen LLC_2_18=rho_2*(1-rho_2^18)/(1-rho_2)
  
 gen LLC_2_4  = rho_2*(1-rho_2^4)/(1-rho_2)
 gen LLC_2_8  = rho_2*(1-rho_2^8)/(1-rho_2)
 gen LLC_2_10 = rho_2*(1-rho_2^10)/(1-rho_2)
 gen LLC_2_12 = rho_2*(1-rho_2^12)/(1-rho_2)
 gen LLC_2_24 = rho_2*(1-rho_2^24)/(1-rho_2)
 
 gen theta_m_r2_18_L   = theta*LLC_5_18/LLC_2_18
 gen theta_m_r2_4_L    = theta*LLC_5_18/LLC_2_4 
 gen theta_m_r2_8_L    = theta*LLC_5_18/LLC_2_8 
 gen theta_m_r2_10_L   = theta*LLC_5_18/LLC_2_10
 gen theta_m_r2_12_L   = theta*LLC_5_18/LLC_2_12
 gen theta_m_r2_24_L   = theta*LLC_5_18/LLC_2_24
 
 gen rho_12=1/(1+0.12)
 gen LLC_12_18=rho_12*(1-rho_12^18)/(1-rho_12)
  
 gen LLC_12_4  = rho_12*(1-rho_12^4)/(1-rho_12)
 gen LLC_12_8  = rho_12*(1-rho_12^8)/(1-rho_12)
 gen LLC_12_10 = rho_12*(1-rho_12^10)/(1-rho_12)
 gen LLC_12_12 = rho_12*(1-rho_12^12)/(1-rho_12)
 gen LLC_12_24 = rho_12*(1-rho_12^24)/(1-rho_12)
 
 gen theta_m_r12_18_L   = theta*LLC_5_18/LLC_12_18
 gen theta_m_r12_4_L    = theta*LLC_5_18/LLC_12_4 
 gen theta_m_r12_8_L    = theta*LLC_5_18/LLC_12_8 
 gen theta_m_r12_10_L   = theta*LLC_5_18/LLC_12_10
 gen theta_m_r12_12_L   = theta*LLC_5_18/LLC_12_12
 gen theta_m_r12_24_L   = theta*LLC_5_18/LLC_12_24

 gen rho_18=1/(1+0.18)
 gen LLC_18_18=rho_18*(1-rho_18^18)/(1-rho_18)
  
 gen LLC_18_4  = rho_18*(1-rho_18^4)/(1-rho_18)
 gen LLC_18_8  = rho_18*(1-rho_18^8)/(1-rho_18)
 gen LLC_18_10 = rho_18*(1-rho_18^10)/(1-rho_18)
 gen LLC_18_12 = rho_18*(1-rho_18^12)/(1-rho_18)
 gen LLC_18_24 = rho_18*(1-rho_18^24)/(1-rho_18)
 
 gen theta_m_r18_18_L   = theta*LLC_5_18/LLC_18_18
 gen theta_m_r18_4_L    = theta*LLC_5_18/LLC_18_4 
 gen theta_m_r18_8_L    = theta*LLC_5_18/LLC_18_8 
 gen theta_m_r18_10_L   = theta*LLC_5_18/LLC_18_10
 gen theta_m_r18_12_L   = theta*LLC_5_18/LLC_18_12
 gen theta_m_r18_24_L   = theta*LLC_5_18/LLC_18_24

 
// Select the robustness check to perfom here:
// ********************************************
 replace theta = theta_m_r18_24_L
// ********************************************
  
 gen theta_4q =. 
 replace theta_4q = 1 if theta < 0
 replace theta_4q = 2 if theta >= 0 & theta < 0.6202
 replace theta_4q = 3 if theta >= 0.6202 & theta < 1.2825
 replace theta_4q = 4 if theta >= 1.2825
 
  collapse(sum) weight*, by(theta_4q)

//forvalues i=1(1)10 {
//forvalues i=1(1)7 { 
forvalues i=1(1)16 { 
 capture gen eta_avg_tmp`i' = eta*weight`i' 
 capture gen theta_avg_tmp`i' = theta*weight`i' 
 capture gen weight_100`i' = round(100*weight`i')
 capture egen eta_avg`i' = sum(eta_avg_tmp`i') 
 capture egen theta_avg`i' = sum(theta_avg_tmp`i')
}  

//drop eta_avg_tmp*
drop theta_avg_tmp*

reshape long weight weight_100 theta_avg, i(theta_4q) j(bs)

keep if weight_100!=0
sort theta_4q
by theta_4q: egen mean_weight_theta=mean(weight_100)
by theta_4q: egen sd_weight_theta=sd(weight_100)
by theta_4q: egen count_weight_theta=count(weight_100)
gen se_weight_theta=sd_weight_theta/(count_weight_theta^0.5)
gen ICl_weight_theta=  mean_weight_theta-1.96* se_weight_theta
gen ICh_weight_theta=  mean_weight_theta+1.96* se_weight_theta

collapse(mean) mean_weight se_weight ICl_weight ICh_weight,by(theta_4q)
	mvencode se_weight,mv(0) over
	replace  ICl_weight = mean_weight if ICl_weight == .
	replace  ICh_weight = mean_weight if ICh_weight == .
	gen tstats = mean_weight/se_weight
	replace tstats = 0 if tstats == .
	order theta mean_weight se_weight ICl_weight ICh_weight tstats 
	sort theta
//outsheet using $pathresults/plot_marginaltheta_4q_2_12_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, comma noquote noname replace 

//Here: extract distribution for each quartile for Excel Table
break
pause


 
//Compute E[m] using heterogeneous r.
//
insheet using $Est_pathresults/Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, clear

 forvalues i=1(1)16 { 
  capture destring weight`i', replace force 
 }

 gen rho_5=1/(1+0.05)
 gen LLC_5_18=rho_5*(1-rho_5^18)/(1-rho_5)
  
 gen rho_2=1/(1+0.02)
 gen LLC_2_18=rho_2*(1-rho_2^18)/(1-rho_2)
  
 gen rho_12=1/(1+0.12)
 gen LLC_12_18=rho_12*(1-rho_12^18)/(1-rho_12)
 
 gen m_2=LLC_2_18/LLC_5_18
 gen m_12=LLC_12_18/LLC_5_18
 gen theta_m_adj = 1
 gen theta_m_2_v    = theta*LLC_5_18/LLC_2_18
 gen theta_m_12_v   = theta*LLC_5_18/LLC_12_18
        
 replace theta_m_adj=theta_m_2_v if theta > m_2
 replace theta_m_adj=theta_m_12_v if theta < m_12		
       

 collapse(sum) weight*, by(theta theta_m_adj)
 
//forvalues i=1(1)10 {
//forvalues i=1(1)7 { 
forvalues i=1(1)16 { 
 capture gen theta_avg_tmp`i' = theta*weight`i' 
 capture egen theta_avg`i' = sum(theta_avg_tmp`i')
 
 capture gen theta_m_adj_avg_tmp`i' = theta_m_adj*weight`i' 
 capture egen theta_m_adj_avg`i' = sum(theta_m_adj_avg_tmp`i')
 
}  

drop theta_avg_tmp*
drop theta_m_adj_avg_tmp*
drop weight*

reshape long theta_avg theta_m_adj_avg, i(theta) j(bs)
	   
collapse(firstnm) theta_avg theta_m_adj_avg,by(bs)
	gen theta_avg_sd = theta_avg
	gen theta_m_adj_avg_sd = theta_m_adj_avg
	gen theta_avg_n = theta_avg
	gen theta_m_adj_avg_n = theta_m_adj_avg
collapse(mean) theta_avg theta_m_adj_avg (sd) theta_avg_sd theta_m_adj_avg_sd  (count)	theta_avg_n  theta_m_adj_avg_n
	gen theta_avg_se = theta_avg_sd / (theta_avg_n^0.5)
	gen theta_m_adj_avg_se = theta_m_adj_avg_sd / (theta_m_adj_avg_n^0.5)
	keep theta_avg theta_m_adj_avg theta_avg_se theta_m_adj_avg_se	   
	   
outsheet using $pathresults/mean_theta_theta_m_adj_Fox_`FOX_model'_mixedFEdemoshrt_WTP_v`subsample'_inc`inc'_sd`set_seed'_cty_tau1_grid_v3_all_bs.csv, comma noquote noname replace 
